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\ The departure of a granular gas in the instable region of parameters from the initial homogeneous 

. cooling state is studied. Results from Molecular Dynamics and from Direct Monte Carlo simulation 

of the Boltzmann equation are compared. It is shown that the Boltzmann equation accurately 
predicts the low density limit of the system. The relevant role played by the parallelization of the 
velocities as time proceeds and the dependence of this effect on the density is analyzed in detail. 
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Q . I. INTRODUCTION 

The study of granular systems has attracted much interest in the last years. The practical importance of these sys- 
tems, which are present in many situations in daily life, together with the richness and complexity of their behavior [1], 
^ •■ have motivated many researches trying to describe and understand the physical mechanisms governing granular flows. 
c/2 ' Of course, the variety of what we call "granular media" makes it necessary to use different theoretical descriptions 
depending on the problem we wish to address. In the context of rapid, low density, granular flows, the application of 
the methods of the kinetic theory for molecular (elastic) systems has proven to be a very useful tool. The starting 
point for this description has been in many cases the extension to inelastic collisions of the Boltzmann equation, which 
^ ■ is derived under the same hypothesis as in elastic systems. From this kinetic equation, closed hydrodynamic equations 
', with explicit expression for the fluxes and the associated transport coefficients (up to second order in the gradients) 
Q ■ for inelastic hard particles have been derived [2,3]. This hydrodynamic picture has been found to provide an accurate 
O ' description of granular systems in very different situations, driven and not driven. The Boltzmann equation has also 
been used to study the velocity distribution of a granular gas, modeled in most of the cases as a system of inelastic 
I hard particles. Of course, the shape of the distribution depends on the state of the system. The simplest possibility is 
^ , the so-called homogeneous cooling state (HCS), the state of a homogeneous, freely evolving granular gas [4]. In that 
' case, the Boltzmann equation is shown to admit a solution whose time dependence can be scaled through its second 
moment. Deviations from gaussianity in the HCS have been quantified by computing the kurtosis of the distribution 

■ [5-7] and by establishing the existence of exponential velocity tails [6,8,9]. The possible solutions of the Boltzmann 
\ equation in the case of a vibrated system in absence of gravity have also been investigated, and the conditions for 
' the existence of a steady solution whose space dependence is scaled out also through the second moment has been 

■ established [10]. 

"p-, , In spite of its extended use and clear success in many cases, the validity of the Boltzmann equation to describe 
granular flows has been questioned since the early developments of the kinetic theory of granular systems. One of the 
^ first objections raised against its use was based on the tendency of these systems to form density clusters. Since it is 
I only valid in the low density limit, the Boltzmann equation is not suitable to describe the cluster evolution. This is 
' O not a fundamental objection, and in fact it can also be raised against the applicability of the Boltzmann equation to 
^ . molecular systems in inhomogeneous states. The value of the density limits indeed the applicability of the Boltzmann 
description, but both in the elastic and the inelastic case. On the other hand, if we start from an homogeneous, 
low density, initial configuration of a granular gas, under conditions that clusters will develop eventually, it can be 
expected that the Boltzmann equation will be valid to describe the first stages of the cluster formation, as long as 
regions with a too large density do not show up. A different objection concerns the validity of the molecular chaos 
; I ■ hypothesis, which is on the basis of the Boltzmann description. As collisions in granular systems tend to make the 
I particle velocities more parallel, velocity correlations may develop from the early stages of the evolution of a granular 
gas, making the molecular chaos hypothesis invalid. This problem was directly addressed by Soto et al [11,12], who 
studied the short-range velocity correlations in the HCS of a granular fluid, concluding that they were not relevant 
in the low density case. Pagonabarraga et al [13] also studied the validity of the molecular chaos hypothesis in a 
homogeneous granular system, but in this case driven by a random force. Again, it was found that, for dilute systems, 
deviations from molecular chaos were not significant. To put this work in a proper context, it should be taken into 
account that the driving force introduced by the authors induced some velocity correlations, as pointed out by them. 
The previous mentioned studies are restricted to homogeneous states of a granular gas. They correspond to 
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very special situations, so the conclusions cannot be extrapolated to more general cases. A different situation was 
investigated recently by Nakanishi [14], who studied the time evolution of the velocity distribution of a freely evolving 
grarndar fluid in conditions such that the HCS is unstable. He found that the kurtosis of the distribution did not 
remain stationary, but evolved towards the Gaussian value from the early stages of the evolution of the system, even 
before the clustering instability shows up. He argues that this is in contradiction with the predictions of kinetic 
theories bascid on the Boltzmann equation, bc;ing a clear indication of the growth of velocity correlations, which 
invalidate the molecular chaos hypothesis. He also fomid that the return to the Gaussian behavior became slower the 
more dilute the system, but as changing the density implies changing the "clustering time" , it was not clear to the 
author how the Boltzmann behavior could be recovered. 

The object of this work is to investigate the possible differences in the behavior of the velocity distribution function 
of a initially dilute, homogeneous, granular gas and the predictions of the Boltzmann equation when the system 
is under conditions such that the HCS is unstable. We will compare the results from Molecular Dynamics (MD) 
simulations at different average densities with those from the Direct Simulation Monte Carlo (DSMC) method, which 
is a method to solve numerically the Boltzmann equation [15]. 

The plan of the paper is as follows: in Section II the kinetic theory predictions for the HCS are briefly discussed, 
while the details of the simulations and the properties to be studied are given in Sec. HI. The results for the DSMC 
and MD simulations are presented in Sees. IV and V. Some final comments and discussion are made in Sec. VI. 



II. KINETIC THEORY FOR THE HOMOGENEOUS COOLING STATE 

Let us consider a a system of N smooth inelastic hard particles, spheres {d = 3) or disks {d = 2), of mass m and 
diameter a. The loss of energy in collisions is characterized by a constant coefficient of normal restitution a, which 
implies that the velocities of two colliding particles i, j before and after the collision are related by 

, 1 + a _ 
Vi = Vi —{^■<t)<t, 

=Vj + ^-(g-cr)o-, (1) 

where the primes denote velocities after the collision, g = v, — Vj is the relative velocity, and a a unit vector joining 
the centers of particles i and j at contact. Let us notice that a < 1, and that the value a = 1 corresponds to elastic 
collisions. The above rule implies that, in each collision, the component of the relative velocity in the direction of a 
is reduced in a factor a: 

g' • ? = -a(g • a) . (2) 

The Boltzmann equation describing the evolution of the velocity distribution /(r, v, t) of a system of freely evolving 
hard particles has the form: 

+ V • /(r,v,i) = JB[r,v|/(r,v,t)] , (3) 

where J7b is the (inelastic) Boltzmann collision operator: 

JB[r,v\f{r,v,t)] = a^-' j dvi j d$e{g-$)ig-$){a-^b-'-l)f{r,v,t)f{r,vut). (4) 

Here is the Heaviside step function and an operator transforming velocities v and vi to its right into their 
precollisional values. As it is the case with elastic particles, the derivation of this equation is based on the molec- 
ular chaos hypothesis, i.e., the factorization of the two-particle distribution function in the precollisional sphere: 
/(^^(xi,X2,f) = /(xi, t)/(x2, i), where x^ = {ri,Vi}. If spatial and/or velocity correlations are present between 
colliding particles, this factorization does not hold and the assimiption of molecular chaos fails. 

When a system of inelastic particles such as the one described evolves freely, its simplest possible state is the HCS. 
It is a homogeneous state with no fluxes, whose temperature T{t), defined as proportional to the average kinetic 
energy, evolves according to Haff 's law [16] , 

m ^ (5) 
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where to is the time characterizing the energy decay. 

When the system is in the HCS, the Boltzmann equation admits a solution /(v, that obeys the scahng law [5-7] 



where Hh is the homogeneous density, and vq the thermal velocity of the system, defined as Vq = y^fceT/m, with 
ks the Boltzmann constant. Therefore, in the HCS all the time dependence of the velocity distribution can be scaled 
out through its second moment. The function (f> is determined from the Boltzmann equation. Although its exact 
expression is not known, it was found that it does not deviate much from a Gaussian [7]. Then, it is sensible to 
expand it using the Sonine polynomials, S^^\ whose explicit expression can be found in [17], 



^(-) = ;^E'^i^^^Hc^), c=^. (7) 

i>o 

Normalization and scaling imply ao = 1 and oi = 0. The coefficient 02 is related to the kurtosis of the distribution 

{v'^) d{d + 2) 
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and its value has been estimated from the Boltzmann equation up to linear order in 02 [5,6]. Both theoretical 
calculations and numerical simulations of the Boltzmann equation [7] show that, for not too inelastic systems, 
is very small. Deviations from gaussianity are important if wc consider the tails of the distribution, which are 
exponential. Again, this has been verified from theoretical arguments based on the Boltzmann equation [8,6] and 
from DSMC simulations [9]. 

The time to characterizing the kinetic energy decay in the HCS has also been computed by using the Boltzmann 
equation [2] , and its expression reads: 

*o =^(")(2T^)rw)(^j 

where C*(q:) is a function that depends only on the coefficient of restitution a, and that, for not too inelastic systems 
reads [2,3] 

C{a)^^^{l-a^). (10) 

In terms of the average number of collisions per particle, r, Half's law takes the form: 

T(t) = T(0)e-2^*^ , (11) 

i.e., the kinetic energy decays exponentially with the number of collisions. 

It is a well known result [18,19] that the homogeneous cooling state of a granular gas is imstablc against long 
wavelength perturbations. In practice, this implies that, if the system size is larger than a critical value, L^, it 
will spontaneously develop spatial inhomogeneities. The value of the critical size has also been determined from the 
Boltzmann equation [3,20], and it reads: 

C,i2^d)nd/2)(,*\''\^^ (12) 



" 27r(<*-3)/2 \^2C* 

Here, \h = 1/ {CdnH(y^'^~^^) is the homogeneous mean free path, Cd = 2\/2 if d = 2 and Cd = 7r-\/2 if d = 3, and 77* 
is a function of a related to the viscosity of a granular gas, whose expression can be found in [3]. Both MD [18,19] 
and DSMC [21] simulations of freely evolving granular gases show that, in the development of inhomogeneities, first 
velocity vortices appear in the system and then the density becomes also very inhomogencous. Of course, once the 
instability has set in, the velocity distribution of the system will no longer be the one of the HCS. Therefore, one 
expects the scaling to fail and 02 to become time dependent. The question is whether the Boltzmann equation is 
valid to describe how the system departs from the HCS. It might happen that, as pointed out by Nakanishi [14], 
the development of velocity correlations are very important from the early stages of this departure, and then the 
Boltzmann equation fails. This is the main point we want to clarify in this work. 
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III. COMPUTER SIMULATIONS 



Given that we want to check the validity of the Boltzmann equation to describe granular gases, we will compare MD 
simulations of initially homogeneous, low density, granular systems with the predictions of the Boltzmann equation. 
Since the analytical solution of the latter is not known, we will use the DSMC method developed by Bird [15] to 
construct numerical solutions of the Boltzmann equation under the desired conditions. It must be remembered that 
this method mimics the dynamics behind the Boltzmann equation by uncouplig the free flow of particles and collisions 
during a small enough time interval. Besides, collisions are always treated as if there were no correlations between 
colliding particles. In practice, this implies that space is divided in cells of size smaller than the mean free path, and 
particles within a cell collide with a probability proportional to their relative velocity. Technical details about the 
DSMC method have been extensively discussed in the literature [15,22]. The only point we want to stress is that, 
when using this method, the simulated system is by definition in the low density limit, where the Boltzmann equation 
is supposed to apply, no matter how many simulated particles there are in a given space region. 

The simulations we will present in the following correspond to a two-dimensional system of freely evolving hard 
disks in a square box of side L, larger than the critical size L,. given by Eq. (12). Periodic boundary conditions will be 
considered in all the cases. In the MD simulations, the event driven algorithm [23] will be used. For our purposes, it is 
important to study the effect of lowering the density on the evolution of the system. For that reason, two different, low 
values, of the density will be studied in the MD simulations, namely nn = 0.05(t^^ and nn = 0.0125cr~^. Nevertheless, 
it must be noticed that, for a given a, changing the density implies changing the critical size of the system, as it 
follows from (12). As a consequence, if the systems have the same size L, the characteristic time governing the growth 
of instabilities will be also changed, and the departure from the HCS will be faster in the denser system. Therefore, 
if we want to have the same "distance to stability" in systems with different densities, their sizes must be changed 
correspondingly, so L/Xh remains the same. The same value of this scaled size should be considered in the DSMC 
simulations, again to expect the same characteristic time governing the departure from the HCS. 

Moreover, we have used a fixed value of the restitution coefficient, a = 0.9, which is not too inelastic, but lies 
outside of what can be considered the quasi-clastic region. For this value of a, Eq. (12) leads to L^. ^ 23.36A/f . As 
we want the system to be well inside the unstable region, we have chosen the system size to be L = 80 A//. Then, for 
riH = 0.05fT~^ the number of particles in the system was N = 16000, while for uh = 0.0125f7~^, A'' = 64000. In the 
DSMC simulations, N = 2.048 • 10^ particles were considered, but it must be reminded that this number has only a 
statistical meaning. In all the cases, we started with the particles homogeneously distributed in the system and with 
a Gaussian velocity distribution. This situation was let to evolve with elastic collisions during several collisions per 
particle, until the system was equilibrated. Then, the inelasticity was switched on, and this time is taken as the origin, 
f = 0, of our simulations. The initial elastic period ensures that the structure of the fluid at t = is the equilibrium 
one. 



A. Studied properties 

As the aim of the paper is to investigate the departure of the velocity distribution of the system from the HCS 
distribution, this will be one of the properties to be studied in the simulations. To be more precise, we will compute 
the second {v"^) and fourth (v**) velocity moments of the total velocity distribution, and from them we will compute 
also 02, given by Eq. (8). As far as the system stays in the HCS, 02 will be constant. It is important to notice that, 
as in Ref. [14], we will study the global velocity distribution of the system, even if inhomogeneities have appeared, 
causing local averages of the physical quantities to be quite different from their global values. 

The growth of inhomogeneities will be controlled by following the evolution of the density n{r,t), and velocity 
u{r,t), fields. To compute them, the system is divided into 30x30 square cells of size Ic ^ 2.7Xh, and properties are 
averaged in each cell. Again, while the system remains in the HCS, n(r, t) should be constant and u(r, t) — (apart 
from statistical noise). Once the instability develops, vortices and also density inhomogeneities will show up. The 
question remains of how the departure of the velocity distribution from the HCS one is related to the development of 
these inhomogeneities. 

The possible failure of the molecular chaos hypothesis will be investigated by the study of several coUisional averages, 
i.e., averages over colliding pairs of particles. These averages contain information about the two-particle distribution 
of colliding particles, which is the function whose factorization is assumed in the molecular chaos hypothesis. The 
first of the coUisional averages we will consider is the pair distribution function at contact, g{a~^). In a clastic system, 
if there are no correlations between colliding particles, as it is in fact assumed by the Boltzmann equation, g{a'^) = 1. 
In the case of a system of inelastic disks, the pair distribution function at contact can be easily computed by using 
[24] (notice that there is a miss-print in the expression for g{cT~^) provided in the cited reference): 
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where we are summing over all collisions 7 taking place in the interval At, and the function implies that we are 
using the pre-coUisional values of the quantities in the sum. When the system is in the HCS, it has been found [24,11] 
that gHCs{(^~^) does not depend on time, and takes the value: 

gHCs{(T+) = ^-^9o{(T+) , (14) 

with go the equilibrium elastic pair distribution function at contact, that in the two dimensional case is accurately 
given by [25] 

, , 7r(25 - 4n7r)n 

^^(-^^ = '+ 4(4 -n.)^ • ^''^ 

In our simulations, time will be discretized so the average over collisions in (13) will be done not over the whole 
simulation time, but over collisions occurring in a given interval. In this way, we can study the evolution of the 
collisional averages and, in particular, of the pair distribution function at contact. Needless to say, g{a~^,t) will only 
be computed in the MD simulations, as in the DSMC method its value is given. 

Another test of the validity of the molecular chaos hypothesis will be provided by the probability distribution of the 
impact parameter, & = § sinx, where x is the angle formed by the impact relative velocity g and the unit vector a. 
In a system of hard disks, and if there are no correlations, the impact parameter is uniformly distributed. Soto et al 
[11] measured this distribution in the HCS and found no deviations from uniformity. Here, we will consider the same 
discretization discussed above to construct the distribution of impact parameters in each time interval, and study its 
possible time dependence. Again, the impact parameter distribution will be studied only in the MD simulations, as 
the DSMC method assumes its uniformity. 

The possible velocity correlations will also be studied. A quantity that has been used to study them is the collisional 
average r{t) defined as 

rW = ;^E (16) 

■'^^ 76 At S I 

with J\fj the number of collisions in the interval At, and i, j the particles involved in collision 7. If there is neither 
macroscopic velocity field nor velocity correlations in the system, it is F = 0. Therefore, in the HCS, non-zero values 
of r are a clear signal of the presence of velocity correlations in the system. When the system leaves the HCS, there 
is another possible reason for these non-zero values. The build-up of velocity vortices implies that u{r, t) is different 
from zero, and then, even if there are no velocity correlations, T will be different from zero. Let us notice that T{t) 
can be mcasiircd both in Molecular Dynamics and DSMC simulations. 

Finally, the distribution of the angle (j) formed by the velocities v^, vj, of colliding particles, will also be computed. 
We are not aware of any analytic expression for this quantity in the elastic case, and, therefore, we will use the 
initial, clastic part, of the simulations to obtain the elastic distribution of the (A angle. The parallclization mechanism 
inherent to inelastic collisions may cause deviations from this behavior, even in the HCS. In any case, this distribution 
shows if there is a predominance of collisions between particles moving more or less parallel in the system. Again, 
this distribution will be measured both in MD and DSMC simulations. 

It is important to remark the different physical information behind the distributions P{b,t) and P{4),t). The 
former is related to the spatial distribution of the incident flux over colliding particles, while the latter contains only 
information about the relative direction of the velocities of colliding particles. 



IV. RESULTS 



The results wc will present correspond, as we have already stated, to MD and DSMC simulations of a freely evolving 
system of hard disks with a = 0.9 and L = 80Ai/, i. e., in conditions such that the HCS is well inside the unstable 
region. In the following, the mass m of the particles will be used as the unit of mass, and the initial kinetic energy as 
the unit of energy. The results we will present have been averaged over several trajectories in all the cases in order to 
improve the statistics. Besides, the time evolution of the system will be expressed in terms of the number of collisions 
per particle, r. As the system is prepared in an initially homogeneous situation, we expect it to stay in the HCS for 



5 



a transient period, until the instability sets in. From, that moment, the different physical properties will depart from 
their HCS values. 

In Fig. 1 we have plotted the evolution of the average kinetic energy per particle, {v'^)/2, in the system. This 
quantity is proportional to the granular temperature as far as there is no macroscopic velocity field, i.e., while u = 0. 
Also plotted is Haff's law, Eq. (11), describing the evolution of this property in the HCS. In all the cases there is an 
initial period in which the energy follows Haff's law and, after that, the average kinetic energy of the system decays 
slower than in the HCS. The departure form Haff's law occurs sooner in the MD simulations than in the DSMC one, 
but it is important to notice that, if we compare the two MD simulations, it occurs sooner in the case of larger density. 
Therefore, it seems that, although there are differences between the behavior of a finite density granular gas and the 
predictions of the Boltzmann equation, they are smaller the smaller the density. Then, it seems sensible to expect 
that the Boltzmann equation predictions provide indeed the correct picture in the analytic low density limit, n — > 0. 

The evolution of 02 is plotted in Figs. 2. In Fig. 2(a), the initial evolution of this quantity is shown: there is an 
initial, very fast, decay from the initial Gaussian, 02 — 0, value to the HCS one, followed by a steady period while 
the velocity distribution of the system remains in the HCS. The complete time evolution of the system is given in 
Fig. 2(b). Again, in all the cases we observe a similar behavior. After the initial steady period, 02 grows, reaches 
a maximum, and afterwards it decays in time. The approximate duration of the steady period, Tst, is Tst ~ 15 for 
Uh = 0.05(7"^, Tst ~ 30 for nn = 0.0125(t~^ and Tst ~ 50 in the DSMC simulation. The position of the maximum, 
r*, is T* ~ 40 for ur = OMa'^, t* ~ 50 for hh = 0.0125cr-2, and r* - 70 in the DSMC simulation. The height 
of the maximum is also larger the larger the density. A first conclusion that can be extracted from Figs. 2 is that 
the qualitative behavior of a2 that follows from the Boltzmann equation is the same found in the MD simulations. 
Besides, the differences between the MD simulations and the DSMC one are smaller the smaller the density in the 
formers, showing again a tendency to the Boltzmann behavior in the limit of very low density. 

Comparison of Figs. 1 and 2 is interesting in order to determine the sensitivity of Haff's law to the exact shape of 
the velocity distribution function. It is found that, for t = Tst, when a2 begins to depart from the steady value, the 
temperature takes the value predicted by Haff's law in all the cases, as it is expected. What it might be surprising is 
that, at the maximum r*, the deviations from the Haff value are relatively small, the temperature being at most 1.3 
times the Haff's value. This is in agreement with simulations of homogeneous systems of inelastic hard disks, that 
show that the evolution of the temperature follows quite approximately Haff's law, unless the fourth moment of the 
velocity distribution is very large as compared with the second one [26] . 

Once the validity of the Boltzmann description as a limit to the behavior of a granular gas for vanishing density has 
been established (at least with regards to the behavior of the second and fourth moment of the velocity distribution) , 
a natural emerging question is which is the mechanism taking the velocity distribution out of the HCS shape. With 
that purpose, the evolution of the density and velocity fields have been studied. It must be noticed that, as we expect 
them to be inhomogencous, and their spatial distribution changes from one realization to another, these fields cannot 
be averaged over different runs. 

Let us consider first the evolution of the density field. In all the cases we observe that the system remains homoge- 
neous during what we have called the steady period, Tgt- Between Tst and r* small density inhomogeneities begin to 
show up, but they are not very significant. The situation changes from r* on: there is a very fast growth of density 
inhomogeneities, and in all the cases, elongated clusters of particles are formed, similar to those observed in previous 
studies of freely evolving granular systems. This behavior of the density field is qu antified in Figs. 3 and 4. In the 
first of the them, the evolution of the dispersion of the density fluctuations, 5 = yj ((n(r, t) — urY) is plotted. Let 
us notice that this quantity is non-zero even in a homogeneous state, due to statistical noise, and its value in the 
homogeneous state depends on the number of particles per cell, which is different in each of the simulations. For that 
reason 5 has been scaled with its value in the elastic part of the simulation, S^i- It is observed in the figures that, in 
all the cases, 5 has not increased much over its elastic value at r = t*, but from then on there is a very sharp growth 
of the density fluctuations as a consequence of the formation of clusters in the system. A similar behavior is exhibited 
by the maximum value of the density, Um, which is plotted in Fig. 4, scaled with the homogeneous density. While 
the system is in a homogeneous state, um is a measure of the statistical noise. Let us remember that the systems 
were divided in 30x30 cells to compute the hydrodynamics fields: as the number of particles is smaller in the denser 
system, the noise in the fields will be larger, and that is the reason why, in the homogeneous part of the evolution, 
nu/nH is larger the larger the density. It must be also pointed out that, once the clustering begins, there is a time 
interval in which the growth of Um can be fitted to an exponential, 

UM ^UH + Ce''^^^-^'^ (17) 

where C and n are constants whose value is not relevant for the present discussion, and s„ ~ 0.13 for n = 0.05ij~^, 
s„ ^ 0.163 for n = 0.0125(7"^, and ,s„ ^ 0.167 in the DSMC simulation. Again, the discrepancies between the 
Boltzmann behavior and the one at finite density are smaller as we lower the density, being the growth rate almost 
the same for the lower density case and the DSMC simulation. 
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In a freely evolving granular system, it has been shown that, when using r as the time variable, the growth rate of 
the scaled transversal velocity mode is 

s±=C- yfc' , (18) 

where fc is a non-dimensional wave- vector defined as A: = 2^ / {^/TrriHCfl), with / the wave-length of the perturbation. 
In a simulation, the maximum allowed wavelength, which is the one that grows faster, \sl = L. For the values of the 
parameters used in this work, the theoretical prediction is s±_ ~ 8.64-10"^. It has been argued [18,27] that the density 
inhomogeneities growth in a freely evolving granular system is a consequence of the non-linear coupling between the 
transversal velocity mode and the other hydrodynamic fields, in particular, the density. If that is indeed the case, the 
growth rate of the density, at least in its first stages, should be 2s_l. Here, 2sj_ ~ 0.173, which is very close to the 
value of Sn found in the DSMC and in the MD lowest density simulation, confirming once again the picture described 
above. This agreement indicates the hydrodynamical character of the density fluctuations in the clustering regime. 

The growth of the scaled velocity field has also been followed in our simulations, and we found that vortices 
begin to develop quite soon in the system. To quantify this, we introduce (5„, the scaled average value of u^, as 
5u = ('^■(r, t)u^(r, t)). Again, due to statistical noise, (5„ will be different from zero when computed from a simulation, 
even if there are no fluxes. For that reason, in Fig. 5 we have plotted 5„ scaled with its value in the elastic part of 
the simulation, 8^ . All the simulations show that the velocity field grows from the very early stages of the evolution. 
In fact, there is a first part of the growth of which can be quite well approximated by an exponential, that 

is common to all the simulations. After that, there is a saturation effect that translates into deviations from the 
exponential growth for larger times, occurring sooner the larger the density. In fact, the saturation occurs for times 
of the order of r* , which is the time when the clustering is triggered. 

V. COLLISIONAL AVERAGES 

The evolution of the coUisional averages in the system has been computed by discretizing the time in intervals 
Ar = 5. Properties arc averaged over collisions taking place in each interval. Also, all the results have been averaged 
over several runs. 

In Fig. 6 we have plotted the evolution of the pair correlation function at contact, g(CT+), from the MD simulations. 
Also included is the theoretical prediction in the HCS for the two values of the density displayed in the figure, 
calculated from Eq. (14). In both cases it is found that, when the velocity distribution begins to depart from the 
HCS one, i.e., at r = Tst, the pair correlation function takes still the HCS value and, in fact, deviations from it are 
quite small even at r = r*. After that, there is a very fast increase of g{(T^) due to the formation of clusters in 
the system. Therefore, the increase in 02 shown in Fig. 2 is not due to the development of positional correlations of 
colliding particles. These correlations do appear, but at rather later times. We have also investigated the evolution 
of the impact parameter distribution, t) and found that in none of the MD simulations discussed in this work it 
deviated significantly from uniformity for the times shown in this manuscript. 

Velocity correlations were investigated through the behavior of F defined in Eq. (16), and of the distribution 
function of the angle formed by the velocities of colliding particles, (j). In Fig. 7, the evolution of F for the three 
simulations is shown. It must be noticed that, in a finite system, even if there are no correlations, F takes a finite 
value that depends on the number of particles, which is different in the three cases. So, even in the elastic part of the 
simulation, F is different in each case, being larger in the denser system, as it has less particles. Besides, the value 
of r is different in a elastic system and in a inelastic one in the same conditions. Therefore, when the inelasticity is 
switched on at t = 0, there is a very fast increase of F to its HCS value. Then, there is a quite steady period that lasts 
longer the lower the density, followed by an almost exponential increase of F. It is found that the growth rate in this 
period depends on the density, becoming larger as n decreases. Finally, there is a slowing down in the increase of F, 
and it seems that it saturates in the end to a value which is larger the larger the density. Also the slowing down begins 
sooner the larger the density. In any case, we find again that by lowering the density in the MD simulations, the 
behavior of the system approaches the one predicted by the Boltzmann equation. One could be tempted to conclude 
that the increase of F is due to the development of correlations between velocities of colliding particles, this happening 
from the early stages of the evolution of the system. But one must be cautious when interpreting this result. Let us 
remember that in our system a velocity field is also being built up from the beginning of the evolution, as it was shown 
in Fig. 5, and this leads to an increase in the value of F that has nothing to do with the failure of the factorization 
of the two-particle distribution function of colliding particles, i.e., with a violation of the molecular chaos hypothesis. 
In fact, the behavior of F displayed in Fig. 7 is quite reminiscent of the one of the fluctuations of the velocity field. 
Also, in the DSMC simulation, it must be taken into account that the molecular chaos hypothesis is assumed in 
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the very basis of the algorithm, so the increase of F in that case cannot be due to the presence of prc-colhsional 
velocity correlations. We conclude then that the growth of F displayed in Fig. 7 is due to the instability of the scaled 
transversal velocity mode, that leads to the formation of vortices and, as a consequence, the velocities of colliding 
particles become more and more parallel. 

The existence of the parallelization mechanism is investigated also by studying the evolution of the distribution 
of the angle formed by the velocities of colliding particles, P{(t>)- In Fig. 8 this distribution is plotted at different 
times for the MD simulation with n = 0.0125cr~^. We have included the distribution at r = 0, which is constructed 
from the elastic part of the simulation. While the system remains in the HCS, the distribution of the (p angle is 
indistinguishable from the clastic one, so there arc no apparent angle correlations in this part of the evolution of 
the system. The situation changes at r ~ Tgt {Tst ~ 30 in this case), when deviations from the elastic distribution 
begin to show up. The relative number of collisions with larger relative angles of the velocity begin to decrease with 
respect to the elastic case. As time proceeds, this tendency becomes stronger, and for the final times considered in 
the simulation, most of the collisions correspond to particles that are moving almost parallel. The distortion of the 
angle distribution begins precisely at the time when the kurtosis of the total velocity distribution begins to depart 
from the HCS value. It seems sensible to conclude then that the behavior of the velocity distribution of the system 
is a consequence of the parallelization mechanism, which is inherent to inelastic collisions, and that, when the system 
is unstable, induces a collective behavior of the velocities of particles, that translates into a departure of the velocity 
distribution from the HCS one. 

The qualitative behavior of P{(p) discussed for n = 0.0125cr~^ also holds for the larger density MD simulation and 
for the DSMC one. Of course, the deviation from the elastic distribution occurs sooner the larger the density, but it 
is also present in the Boltzmann description. To have a clear picture of the parallelization of the velocities of colliding 
particles in the three cases, in Fig. 9 we have plotted the evolution of {(f)) in the three simulations. The value of 
this quantity in an equilibrium, elastic fluid, is {(f>) ~ 0.597r(~ 107°). In the three cases, the system remains with 
the elastic value of {(j)) during what we called the steady period, and then it begins to decrease with time, indicating 
that there is a tendency to have more collisions with velocities that form a small angle. It must be also remembered 
that, in the MD simulations, the impact parameter distribution did not show deviations from uniformity for the times 
considered in this work. 



VI. DISCUSSION 

The simulation of a system on freely evolving hard disks in conditions such that the HCS is unstable shows that 
the Boltzmann description provides a valid picture for the behavior of a low density granular gas. This has been 
established by showing that there arc no qualitative differences between the results obtained in low density MD 
simulations and those that follow by using the DSMC method. It has been shown that, when the density is lowered 
in the suitable way in the MD simulations, i.e., leaving the characteristic time for the development of instabilities 
unchanged, the behavior of the system tends to the Boltzmann one. Nevertheless, it is also true that the deviations 
from the Boltzmann behavior are larger in a system in these conditions of instability than in a stable one, for the 
same values of the density and restitution coefficient. The reason for this seems to be that, when the HCS in unstable, 
there is a parallelization of velocities of colliding particles that is more efficient the higher the density, although it is 
also present in the Boltzmann description. 

The role of the spatial correlations between colliding particles has been investigated in this work by studying the 
pair distribution at contact and the probability distribution of impact parameters. The former only deviate from 
the HCS value when density inhomogeneities are already developed in the system, while the latter remains always 
uniform. This implies that, in a low density granular gas, the development of spatial correlations does not play a 
significant role in the early departure from the HCS. 

The above results seem to indicate that, when trying to extend the inelastic Boltzmann equation to finite higher 
density, the effect of velocity correlations between colliding particles must be incorporated. At least, in the physical 
situation considered in this paper, velocity correlations become important quite before the system develops significant 
spatial correlations, as measured by the pair distribution function at contact. It is sensible to expect that this effect 
increases as the inelasticity of the system increases. If this picture were right, kinetic equations for inc^lastic dense gases 
should not be based on Enskog-like equations, taking into account only spatial correlations, but new approximations 
incorporating the effect of velocity correlations in the precollisional sphere are needed. This implies a rather strong 
departure from the traditional methods of kinetic theory for elastic systems. 
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FIG. 1. Evolution of the average kinetic energy for the simulations discussed in the text. The units are chosen such that the 
initial value is equal to one. Also shown is the theoretical prediction for the HCS, i.e., Haflf's law. Time is measured in average 
number of collisions per particle. 

FIG. 2. Time evolution of the dimensionless coefficient a2 for the simulations discussed in the main text. The short time 
behavior, (a), and the complete evolution, (b), are displayed in diflferent plots for the sake of clarity. 

FIG. 3. Evolution of the average value of the density fluctuations, scaled with their value in the elastic case, for the simulations 
discussed in the main text. 

FIG. 4. Evolution of the maximum value of the scaled density for the simulations discussed in the main text. The symbols 
are from the simulations, and the dotted line is just a guide to the eye. 

FIG. 5. Evolution of the fluctuations of the velocity field, scaled with its value in the elastic case, for the simulations discussed 
in the main text. 
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FIG. 6. Evolution of the pair distribution function at contact, g{cr~^), obtained in the two MD simulations discussed in the 
main text. The horizontal lines are the theoretical prediction for this function in the HCS at n = 0.0125cr~^ and n = 0.05cr~^, 
from bottom to top. 



FIG. 7. Evolution of F for the simulations discussed in the main text. 

FIG. 8. Probability distribution of the angles of the velocities between colliding particles, P(0), from the MD simulation 
with n = 0.0125cT"^ . 

FIG. 9. Evolution of the average value of the angle between the velocities of colliding particles, (0), for the simulations 
discussed in the main text. 
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